* This file creates multiple imputed datasets 

global name "mult_impu"
global main "C:/Users/silvio/Documents/CVR/ryp"
global output "$main/output"  
global data "${main}/data"  
global odata "${data}/Intermuestra1-2"  

global code "${main}/code"  
global coord "$data/coordinates"  
global geocoord "$data/geocoord"  

cd $main
capture log close
log using "${output}/s1_${name}.log", replace
clear

set more off
use "${data}/cvrextrapo", clear
set matsize 1392

scalar nstrata=58
set seed 1234

matrix y=J(nstrata*3*8,4,0)
scalar ij=1

local i=1
while `i'<=nstrata{

scalar perpeh=1
while perpeh<=3{
local j=1
while `j'<=8{
matrix y[ij,1]=ij
matrix y[ij,2]=perpeh
matrix y[ij,3]=`j'
matrix y[ij,4]=`i'
scalar ij=ij+1
local j=`j'+1
}
scalar perpeh=perpeh+1
}
local i=`i'+1
}

matrix ymix=y
matrix mfour=J(nstrata,1,-4)




gen miss1=1-miss
gen ones=1

mlogit perpe2 sexon situa anohechos age cvr odh dp ones, noconst base(3)
mfx compute, predict(outcome(1))
mfx compute, predict(outcome(2))

matrix b=e(b)
matrix V=e(V)
matrix b=b[1,1..16]
matrix V=V[1..16,1..16]
matrix Z=J(1,16,0)

matrix C=cholesky(V)

mlogit perpe2 sexon situa anohechos cvr odh dp ones, noconst base(3)
mfx compute, predict(outcome(1))
mfx compute, predict(outcome(2))

matrix ba=e(b)
matrix Va=e(V)
matrix ba=ba[1,1..14]
matrix Va=Va[1..14,1..14]
matrix Za=J(1,14,0)

matrix Ca=cholesky(Va)

matrix uno=J(1,1,0)
matrix tres=J(3,1,0)

tab perpe2 if j~=.,matcell(k1)
tab perpe2 source if j~=., matcell(k2)
matrix k1=uno,k1'
matrix ktot=k1
matrix k2=tres,k2
matrix ktot2=k2



scalar mit=1
while mit<=20{



gen exp1=exp(sexon*b[1,1]+situa*b[1,2]+anohechos*b[1,3]+age*b[1,4]+cvr*b[1,5]+odh*b[1,6]+ dp*b[1,7]+b[1,8])
 gen exp2=exp(sexon*b[1,9]+situa*b[1,10]+anohechos*b[1,11]+age*b[1,12]+cvr*b[1,13]+odh*b[1,14]+ dp*b[1,15]+b[1,16])
 gen pr1=exp1/(1+exp1+exp2)
 gen pr2=exp2/(1+exp1+exp2)
 gen pr3=1-pr1-pr2
gen pr12=pr1+pr2
gen perpred=3
gen random=runiform()
replace perpred=1 if random<=pr1
replace perpred=2 if random>pr1 & random<=pr12
replace perpred=perpe2 if miss1==1
drop exp1 exp2 pr1 pr2 pr3 pr12 random


gen exp1=exp(sexon*ba[1,1]+situa*ba[1,2]+anohechos*ba[1,3]+cvr*ba[1,4]+odh*ba[1,5]+ dp*ba[1,6]+ba[1,7])
 gen exp2=exp(sexon*ba[1,8]+situa*ba[1,9]+anohechos*ba[1,10]+cvr*ba[1,11]+odh*ba[1,12]+ dp*ba[1,13]+ba[1,14])
 gen pr1=exp1/(1+exp1+exp2)
 gen pr2=exp2/(1+exp1+exp2)
 gen pr3=1-pr1-pr2
gen pr12=pr1+pr2
gen perpreda=3
gen random=runiform()
replace perpreda=1 if random<=pr1
replace perpreda=2 if random>pr1 & random<=pr12
replace perpreda=perpe2 if miss1==1
replace perpred=perpreda if perpreda~=. & perpred==.

drop perpreda

local i=1
while `i'<=16 {
matrix Z[1,`i']=invnorm(uniform())
local i=`i'+1
}
local i=1
while `i'<=14 {
matrix Za[1,`i']=invnorm(uniform())
local i=`i'+1
}

*matrix u=matuniform(1,16)
matrix b=b+Z*C'
matrix ba=ba+Za*Ca'

scalar perpeh=1
quietly: tab j source if perpred==perpeh, nol matcell(w) matrow(z)
run "$code\reconv"
matrix x1=x

scalar perpeh=2
quietly: tab j source if perpred==perpeh, nol matcell(w) matrow(z)
 run "$code\reconv"
matrix x2=x

scalar perpeh=3
quietly: tab j source if perpred==perpeh, nol matcell(w) matrow(z)
 run "$code\reconv"
matrix x3=x




*log on
di mit
tab perpred if j~=.,matcell(k1)
tab perpred source if j~=., matcell(k2)
*log off
matrix uno[1,1]=mit
matrix tres[1,1]=mit
matrix tres[2,1]=mit
matrix tres[3,1]=mit
matrix k1=uno,k1'
matrix ktot=ktot \ k1
matrix k2=tres,k2
matrix ktot2=ktot2 \ k2



matrix dymi=x1,mfour,x2,mfour,x3,mfour
matrix dymi=vec(dymi')

if mit==1{
matrix ymi=dymi
}
else {
matrix ymi=ymi,dymi
}

drop exp1 exp2 pr1 pr2 pr3 pr12 perpred random 

di mit
scalar mit=mit+1
}

matrix ymi=ymi,ymix
clear
svmat ymi, names(y)
ren y21 n
ren y22 perpe
ren y23 k
ren y24 i

order n perpe k i y1-y20
save "${data}/${name}", replace
clear
capture log close
